#if !defined(TYPELIB_OPERATIONS_H)
#define TYPELIB_OPERATIONS_H

#include <math.h>

//!!!Operations should NOT be in the namespace 'types',
//!!!otherwise all their specializations have to be defined in this namespace also,
//!!!and this is inappropriate.
//namespace types {

template <class C> 
inline double scalar_product(const C & x, const C & y)
	{ return x * y; } //if you want something else then create specialization

template <class C>
inline double square(const C & x)
	{ return scalar_product(x, x); }

template <class C> 
inline double abs(const C & x)
	{ return sqrt(square(x)); }

template <>
inline double abs(const int & x)
	{ return x < 0 ? - x : x; }

template <>
inline double abs(const float & x)
	{ return fabs(x); }

template <>
inline double abs(const double & x)
	{ return fabs(x); }

template <class T> 
inline void inverse(T & t)
	{ t = -t; }

///2-point linear interpolation (on line)
template <class T, class F> 
inline T linear(const T & x0, const T & x1, F a)
	{ return x0 * (1 - a) + x1 * a; }

///2-point linear interpolation (on line); inplace version
template <class T> 
inline void linear2(const T & x0, const T & x1, double a, T & x)
	{ x = x0 * (1 - a) + x1 * a; }

///3-point linear interpolation (on plane)
template <class T, class F> 
inline T linear(const T & x0, const T & x1, const T & x2, F a1, F a2)
	{ return x0 * (1 - a1 - a2) + x1 * a1 + x2 * a2; }

///returns weight in linear interpolation of '0'-node for the 'p'-point in an triangle
template <class T>
inline double linear_trinode_weight(const T & x0, const T & x1, const T & x2, const T & p)
	{ return outer_product(x2-x1, p-x1) / outer_product(x2-x1, x0-x1); }

template <class T, class F> 
inline T bilinear(const T & x00, const T & x01, const T & x10, const T & x11, F a0, F a1) {
	return	x00 * (1 - a0) * (1 - a1) +
			x01 * (1 - a0) * a1 +
			x10 * a0 * (1 - a1) +
			x11 * a0 * a1;
}

//} // types

#endif // TYPELIB_OPERATIONS_H
